///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file determinants.cpp
///
/// \brief   implementation for functions associated with determinants
///
///////////////////////////////////////////////////////////////////////////////////////////////////

// #include <src/pcps.h>
// #include <src/determinants.h>
// #include <src/mantissa_exponent.h>
// #include <src/blas_lapack.h>
// #include <src/matrix_inverse_lu.h>
// #include <src/random.h>

// void get_nonzero_in_u(const int n, int & count, int & pos, const double * const u) {
//   count = 0;
//   pos = -1;
//   for (int i = 0; i < n; i++) {
//     const int q = ( u[i] != 0.0 );
//     count += q;
//     pos += q*(i-pos);
//   }
// }

// void get_nonzero_in_u(const int n, int & count, int & pos, const std::complex<double> * const u) {
//   count = 0;
//   pos = -1;
//   for (int i = 0; i < n; i++) {
//     const double x = pcps::square_norm(u[i]);
//     const int q = ( x != 0.0 );
//     count += q;
//     pos += q*(i-pos);
//   }
// }

// ///////////////////////////////////////////////////////////////////////////////////////////////////
// /// \brief   computes the ratio of the determinant of an invertable matrix A to the determinant
// ///          of A plus a rank one update
// ///
// ///    ratio = det( A + u v^T ) / det(A)
// ///
// ///    Where A is nxn and u and v are nx1.
// ///
// /// \param[in]     n        the dimension of the matrix A
// /// \param[in]     A_inv    the inverse of the invertable matrix A, size n*n
// /// \param[in]     u        size n
// /// \param[in]     v        size n
// /// \param[out]    w        workspace, size n
// ///
// /// \return the determinant ratio
// ///
// ///////////////////////////////////////////////////////////////////////////////////////////////////
// template <class T> T pcps::determinant_ratio_rank_1(const int n,
//                                                     const T * const A_inv,
//                                                     const T * const u,
//                                                     const T * const v,
//                                                     T * const w)
// {
// 
//   int u_nonzero, v_nonzero, x;
// 
//   // check how many non-zero elements u has and get the index of a non-zero element
//   get_nonzero_in_u(n, u_nonzero, x, u);
// 
//   // ratio is one if u is all zeros
//   if ( u_nonzero == 0 )
//     return pcps::unity<T>();
// 
//   // simple calculation if u has only one nonzero element
//   if ( u_nonzero == 1 ) {
//     //std::cout << "used u as elementry vector ";
//     return pcps::unity<T>() + u[x] * pcps::xdot(n, v, 1, A_inv + x, n);
//   }
// 
//   // check how many non-zero elements v has and get the index of a non-zero element
//   get_nonzero_in_u(n, v_nonzero, x, v);
// 
//   // ratio is one if v is all zeros
//   if ( v_nonzero == 0 )
//     return pcps::unity<T>();
// 
//   // simple calculation if v has only one nonzero element
//   if ( v_nonzero == 1 ) {
//     //std::cout << "used v as elementry vector ";
//     return pcps::unity<T>() + v[x] * pcps::xdot(n, u, 1, A_inv + x*n, 1);
//   }
// 
//   // in the general case, we must compute a matrix-vector product
//   pcps::xgemm('T', 'N', n, 1, n, pcps::unity<T>(), A_inv, n, u, n, pcps::zero<T>(), w, n);
// 
//   // return the ratio
//   return pcps::unity<T>() + pcps::xdot(n, v, 1, w, 1);
// 
// }
// 
// template double pcps::determinant_ratio_rank_1(const int n,
//                                                const double * const A_inv,
//                                                const double * const u,
//                                                const double * const v,
//                                                double * const w);
// template std::complex<double> pcps::determinant_ratio_rank_1(const int n,
//                                                              const std::complex<double> * const A_inv,
//                                                              const std::complex<double> * const u,
//                                                              const std::complex<double> * const v,
//                                                              std::complex<double> * const w);

// template<class T> void pcps::test_determinants(const int n) {
// 
//   std::vector<T> S(n*n, pcps::zero<T>());
//   for (int i = 0; i < n; i++)
//   for (int j = 0; j < n; j++) {
//     S[i*n+j] = 0.001 * pcps::random_number<T>() - 0.0005*(pcps::unity<T>() + pcps::imaginary_unity<T>());
//   }
//   for (int i = 0; i < n; i++)
//     S[i*n+i] = pcps::unity<T>() * ( 0.8 + 0.4 * pcps::random_number<double>() );
// 
//   std::vector<T> S_inv(S);
//   std::vector<T> S_save(S);
// 
//   std::vector<T> twork(S);
//   std::vector<int> iwork0(n);
//   std::vector<int> iwork1(n);
// 
//   pcps::MantExp<T> d = pcps::matrix_determinant_lu(n, &S.at(0), &iwork0.at(0));
// 
//   T waste;
//   pcps::matrix_inverse_lu(n, waste, &S_inv.at(0), &twork.at(0), &iwork0.at(0), &iwork1.at(0));
// 
//   std::cout << "initial determinant:" << std::endl;
//   std::cout << pcps::format_number("%20.12f", d.get()) << std::endl;
//   std::cout << std::endl;
// 
//   // general u and v
//   {
// 
//     S = S_save;
// 
//     std::vector<T> u(n);
//     for (int i = 0; i < n; i++)
//       u[i] = pcps::unity<T>() * ( -0.5 + 0.5 * pcps::random_number<double>() );
// 
//     std::vector<T> v(n);
//     for (int i = 0; i < n; i++)
//       v[i] = pcps::unity<T>() * ( -0.5 + 0.5 * pcps::random_number<double>() );
// 
//     std::cout << "quick ratio:  ";
//     std::cout << pcps::format_number("%20.12f", pcps::determinant_ratio_rank_1(n, &S_inv.at(0), &u.at(0), &v.at(0), &twork.at(0)) )
//               << std::endl;
// 
//     for (int i = 0; i < n; i++)
//     for (int j = 0; j < n; j++)
//       S[i*n+j] += u[i] * v[j];
// 
//     std::cout << "slow ratio:   ";
//     std::cout << pcps::format_number("%20.12f", ( pcps::matrix_determinant_lu(n, &S.at(0), &iwork0.at(0)) / d ).get())
//               << std::endl;
// 
//     std::cout << std::endl;
// 
//   }
// 
//   // row change
//   {
// 
//     S = S_save;
// 
//     std::vector<T> u(n, pcps::zero<T>());
//     u[pcps::global_lcg() % n] = pcps::unity<T>() * ( 0.5 + 0.5 * pcps::random_number<double>() );
// 
//     std::vector<T> v(n);
//     for (int i = 0; i < n; i++)
//       v[i] = pcps::unity<T>() * ( -0.5 + 0.5 * pcps::random_number<double>() );
// 
//     std::cout << "quick ratio:  ";
//     std::cout << pcps::format_number("%20.12f", pcps::determinant_ratio_rank_1(n, &S_inv.at(0), &u.at(0), &v.at(0), &twork.at(0)) )
//               << std::endl;
// 
//     for (int i = 0; i < n; i++)
//     for (int j = 0; j < n; j++)
//       S[i*n+j] += u[i] * v[j];
// 
//     std::cout << "slow ratio:   ";
//     std::cout << pcps::format_number("%20.12f", ( pcps::matrix_determinant_lu(n, &S.at(0), &iwork0.at(0)) / d ).get())
//               << std::endl;
// 
//     std::cout << std::endl;
// 
//   }
// 
//   // column change
//   {
// 
//     S = S_save;
// 
//     std::vector<T> u(n);
//     for (int i = 0; i < n; i++)
//       u[i] = pcps::unity<T>() * ( -0.5 + 0.5 * pcps::random_number<double>() );
// 
//     std::vector<T> v(n, pcps::zero<T>());
//     v[pcps::global_lcg() % n] = pcps::unity<T>() * ( 0.5 + 0.5 * pcps::random_number<double>() );
// 
//     std::cout << "quick ratio:  ";
//     std::cout << pcps::format_number("%20.12f", pcps::determinant_ratio_rank_1(n, &S_inv.at(0), &u.at(0), &v.at(0), &twork.at(0)) )
//               << std::endl;
// 
//     for (int i = 0; i < n; i++)
//     for (int j = 0; j < n; j++)
//       S[i*n+j] += u[i] * v[j];
// 
//     std::cout << "slow ratio:   ";
//     std::cout << pcps::format_number("%20.12f", ( pcps::matrix_determinant_lu(n, &S.at(0), &iwork0.at(0)) / d ).get())
//               << std::endl;
// 
//     std::cout << std::endl;
// 
//   }
// 
// }
// 
// template void pcps::test_determinants<double>(const int);
// template void pcps::test_determinants<std::complex<double> >(const int);
